/*
    Copyright 2010-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/** \file 

    Declaration of the unstructured mesh grid classes used to
    interface with Arepo. NOTE that because this file pulls in the
    Arepo declarations, which contain many DEFINE's, this file should
    be included AFTER other files, because otherwise the risk of macro
    clashes is large. */

// $Id$

#ifndef __arepo_grid__
#define __arepo_grid__

namespace mcrx {
  template <typename> class arepo_grid;
  template <typename> class arepo_grid_cell;
};

#include "config.h"

#ifdef WITH_AREPO

#include "mcrx-types.h"
#include "ray.h"
#include <vector>
#include "mcrx-debug.h"
#include <cmath>
#include "blitz/numinquire.h"
#include "mcrx-arepo.h"
#include "mcrx-debug.h"
#include "mono_poly_abstract.h"

// KLUDGE: we don't need these two here but they are included by arepo
// and they won't compile in C-linkage mode
#include "mpi.h"
#include "gmp.h"
extern "C" {
#include "proto.h"
#include "voronoi.h"

}
namespace arepo {
  extern tessellation* T;
};

/** Class representing an unstructured mesh grid generated by
    Arepo. This is only an interface class, the grid is actually
    handled by Arepo itself by using its datastructures. Note that
    this means there can not be more than one arepo_grid object, as
    the Arepo data structures are global. */
template <typename cell_data_type>
class mcrx::arepo_grid {
public:
  typedef mcrx::T_float T_float;
  typedef arepo_grid T_grid;
  typedef cell_data_type T_data;
  typedef arepo_grid_cell<T_data> T_cell;
  typedef typename T_cell::cell_tracker T_cell_tracker;
  friend  class arepo_grid_cell<T_data>;
  friend class T_cell::cell_tracker;

protected:
  /// The cell objects are kept in this vector
  std::vector<T_cell> cells_;

  /// The Arepo point array (positions of mesh points).
  point* DP;
  /// The Arepo tetrahedron array.
  tetra* DT;
  /// The length of the Arepo point array.
  int Ndp;

public:
  arepo_grid();

  virtual ~arepo_grid() {};

  T_cell_tracker locate (const vec3d &p, bool=false);

  std::pair<T_cell_tracker, T_float> 
  intersection_from_without (const ray_base& ray);

  /// Returns the number of cells in the grid.
  int n_cells () const {return N_gas;};
  /// Returns cell i (this is the same as the cell with SphP index i).
  T_cell& cell(int i) {assert(i>=0); assert(i<cells_.size()); return cells_[i];};
  /// Returns cell i (this is the same as the cell with DP index i).
  const T_cell& cell(int i) const {return cells_[i];};

  /** Returns the number in the grid enumeration (the DP index) of the
      specified cell. */
  int get_cell_number(const T_cell_tracker& c) const {
    assert(c); return c.cell_->dp_idx_;};

  // This is a dummy because it doesn't make sense for the Arepo grid.
  std::vector<bool> get_structure () const {assert(0);};

  typedef typename std::vector<T_cell>::iterator iterator;
  typedef typename std::vector<T_cell>::const_iterator const_iterator;

  const_iterator begin () const { return cells_.begin(); };
  iterator begin () { return cells_.begin(); };
  const_iterator end () const { return cells_.end(); };
  iterator end () { return cells_.end(); };
};


/** Class representing a Voronoi cell in the Arepo mesh. The actual
    mesh info is stored by Arepo itself, so this just interfaces with
    that. In this context, a "cell" is a Voronoi polygon and uniquely
    defined by its mesh point. */
template <typename cell_data_type> class mcrx::arepo_grid_cell {
public:
  typedef mcrx::T_float T_float;
  typedef arepo_grid_cell<cell_data_type> T_cell;
  typedef cell_data_type T_data;
  typedef arepo_grid<T_data> T_grid;
  friend  class arepo_grid<T_data>;
  friend int main(int, char**);

protected:
  /// The grid this cell is a member of.
  T_grid* g_;

  /** The index into the Arepo DP array of the primary cell that
      corresponds to this mesh cell. Negative corresponds to an
      invalid cell. */
  int dp_idx_;

  /// Pointer to the cell data object.
  T_data* data_;

  arepo_grid_cell(T_grid& g, int p) : g_(&g), dp_idx_(p), data_(0) {};
public:
  arepo_grid_cell() : g_(0), dp_idx_(-1), data_(0) {};

  operator bool() const {return dp_idx_>=0;};

  /// Returns the pointer to the cell data object.
  T_data* data() const {assert (data_); return data_;};
  /** Sets the celldata pointer, deleting the object previously
      pointed to first. */
  void set_data(T_data* d) {if (data_) delete data_; assert (d); data_ = d;};

  /** Sets the celldata pointer, without deleting the object
      previously pointed to. This is used when the current celldata
      was created using placement (this is not optimal). */
  void set_data_nodelete(T_data* d) {data_ = d;};

  const T_grid* host_grid () const {assert (g_); return g_;};
  T_grid* host_grid () {assert (g_); return g_;};

  /** Returns the index into the cells_ array for this cell (which is
      also the index into the SphP Arepo array.) */
  ptrdiff_t index() const {const T_cell* c=&g_->cells_[0]; return this-c;};

  /// Returns the index of the primary cell in the Arepo DP array for this cell.
  int primary_cell() const { return arepo::primary_cells[dp_idx_];};

  /** Returns the position of the primary mesh point associated with
      this cell, converted to Sunrise coordinates. */
  vec3d position() const { return arepo::get_point(primary_cell()); };

  class cell_tracker;
  friend class cell_tracker;
  typedef cell_tracker T_cell_tracker;

  /// Returns the volume of the cell in Sunrise units.
  T_float volume() const {
    using namespace arepo;
    const int sphp_idx= sph_index(dp_idx_);
    return SphP[sphp_idx].Volume*lcon*lcon*lcon; };
};
  

/** This class abstracts the concept of a "current cell" for the
    radiation transfer code. This class knows how to step through a
    Voronoi cell. This is the base class for the cell_tracker and
    const_cell_tracker classes. */
template<typename T_data>
class mcrx::arepo_grid_cell<T_data>::cell_tracker {
public:
  typedef arepo_grid_cell<T_data> T_cell;
  typedef typename T_cell::T_cell_tracker T_cell_tracker;
  friend class arepo_grid<T_data>;
  typedef arepo_grid_cell<T_data>* T_ptr;

private:

  /// The DP index of the Voronoi cell we are within
  int dp_idx_;
  /// The DP index of the cell we came from.
  int previous_cell_;

  /** The pointer to the primary cell we are within. Because of
      periodic BC, several DP cells point to the same primary cell. */
  arepo_grid_cell<T_data>* cell_; 

  int dp_index() const { return dp_idx_; };

  T_cell* cell() const { return cell_; };

public:

  cell_tracker() : 
    dp_idx_(-1), previous_cell_(-1), cell_(0) {};

  cell_tracker(int dp, T_ptr c) : 
    dp_idx_(dp), previous_cell_(-1), cell_(c) {};

  cell_tracker(const cell_tracker& rhs) :
    dp_idx_(rhs.dp_idx_), previous_cell_(rhs.previous_cell_), 
    cell_(rhs.cell_) {};

  ~cell_tracker() {};
  
  /** Assignment operator copies everything. */
  cell_tracker& operator=(const cell_tracker& rhs) {
    dp_idx_ = rhs.dp_idx_;
    previous_cell_ = rhs.previous_cell_;
    cell_ = rhs.cell_;
    return *this;
  }

  /** Equality comparison must use the dp_idx, since several Voronoi
      cells point to the same cell object. */
  bool operator==(const cell_tracker& rhs) const {
    return (dp_idx_ == rhs.dp_idx_); };

  /// Converts to a bool indicating whether cell is valid or not.
  operator bool() const {return cell_ != 0;};

  /// Resets the tracker to be invalid.
  void reset() {dp_idx_ = -1; cell_ = 0; previous_cell_ = -1; };

  /** Resets the cached intersection calculation. This MUST be called
      if the intersection routines are called repeatedly and the ray
      has changed direction or has been moved transverse to its
      direction. */
  void reset_intersections() {
    previous_cell_=-1; };

  T_float 
  intersection_from_within (const ray_base& r, T_float max_len, T_densities& dn);

  /** Converts a fractional column density (optical depth) through the
      cell to a fractional physical distance through the cell. This is
      used when determining where in the cell a scattering takes
      place. */
  T_float column_to_distance_fraction(T_float fn, const vec3d& pos, 
				      const vec3d& dir, T_float L) const {
    return arepo::column_to_distance_fraction(dp_idx_, 
					      fn, pos, dir, L); };

  typename T_cell::T_data* data() const {
    return cell_->data(); };

  /// Checks whether a position is within the cell. 
  bool contains(const vec3d& p) const;

  /// Asserts that a position is within the cell. Used for testing.
  void assert_contains(const vec3d& p) const {
    arepo::assert_contains(dp_idx_, p); };

};


/** Constructor creates cell vector. Before this is called, Arepo must
    have been called to load the snapshot. We don't want to load the
    snapshot here because we might want more than one grid referring
    to the Arepo data. */
template <typename cell_data_type>
mcrx::arepo_grid<cell_data_type>::arepo_grid ()
{
  DP = arepo::T->DP;
  Ndp = arepo::T->Ndp;
  
  // The Ndp mesh cells only refer to N_gas actual cells with hydro
  // quantities. Those are the ones we want.
  cells_.resize(N_gas);
  for(int i=0; i < Ndp; ++i)
    if(arepo::primary_cell(i)) {
      assert(DP[i].task==ThisTask);
      cells_[DP[i].index] = T_cell(*this, i);
    }
}


/** Returns the arepo_grid_cell object that contains the specified
    point. It uses the arepo::locate function. Note that the
    accept_outside parameter doesn't do anything here. \todo Is this
    bad? */
template <typename cell_data_type>
typename mcrx::arepo_grid<cell_data_type>::T_cell_tracker
mcrx::arepo_grid<cell_data_type>::locate (const vec3d & p, bool)
{
  // most of the work is done by this function which has no knowledge
  // of anything outside of Arepo.
  const int dp_idx = arepo::locate(p);

  if(dp_idx>=0) {
    // found 
    const int sphp_idx = arepo::sph_index(dp_idx);
    return T_cell_tracker(dp_idx, &cells_[sphp_idx]);
  }
  else
    // no cell found
    return T_cell_tracker();
}


/** Checks if the ray hits the Arepo box and finds the intersection
 point.  Returns a cell_tracker object representing the cell the
 ray will enter if it hits the box as well as the distance between the
 ray position and the point where the ray enters the box.  If the ray
 misses, NaN and an invalid tracker_info is returned. If the ray is
 already inside the box, the cell containing the ray is returned. */
template <typename cell_data_type>
std::pair<typename mcrx::arepo_grid<cell_data_type>::T_cell_tracker, mcrx::T_float>
mcrx::arepo_grid<cell_data_type>::intersection_from_without (const ray_base& ray)
{
  using arepo::lcon;

  // First test whether the ray already is inside the Arepo box.
  T_float l=0;
  ray_base temp_ray(ray);

  if(!arepo::in_arepo_box(ray.position())) {
    // It is not. Test if it will hit the box.
    l = ray.intersect_cube(arepo::arepomin, arepo::arepomax, false);
    if(l!=l) {
      DEBUG(2,std::cout << "Intersection from without: missed Arepo box" << std::endl;);
      // we missed the box
      return make_pair(T_cell_tracker(), l);
    }
    // It will hit the box. Propagate the temporary ray to the point on the face.
    temp_ray.propagate(l);
  }

  // Find the cell the ray will enter / is in.
  T_cell_tracker t = locate(temp_ray.position(), true);
  if(!t) {
    DEBUG(1,std::cerr << "Error: intersection_from_without got no valid cell even though ray is in box!" << std::endl;);
    return make_pair(t, blitz::quiet_NaN(T_float()));
  }
    
  DEBUG(2,std::cout << "Intersection from without: hit valid cell " << t.dp_index() << " after " << l << std::endl;);
  return make_pair(t, l);
}



/** Returns the result of propagating the ray through the Voronoi cell
    and into the next, subject to a maximum propagation distance of
    max_len. The propagated length is returned and dn is set to the
    column densities in the propagated segment. The ray must be INSIDE
    this cell for this to be called. */
template<typename T_data>
//template<typename T_ptr>
mcrx::T_float
mcrx::arepo_grid_cell<T_data>::cell_tracker::
intersection_from_within (const ray_base& r, T_float max_len, T_densities& dn)
{
  using namespace arepo;

  assert(cell_);
  arepo::assert_contains(dp_index(), r.position());

  point *DP = cell_->host_grid()->DP;

  const int sphp_idx = sph_index(dp_index());
  DEBUG(2,std::cout << "Intersection from within at pos " << r.position() << ", cell " << dp_index() << " sph cell " << sphp_idx << std::endl;);

  // To find the intersection with the next cell, we use
  // find_next_cell which deals with the Arepo data and returns the DP
  // index of the next cell and the distance to the face.
 
  T_float len;
  const int next = arepo::find_next_cell(r, dp_index(), previous_cell_, len);

  // Get distance from current position to the next intersection point
  // (or to the max length point if we would pass it)
  const T_float dl = std::min(max_len - r.length(), len);
  assert(dl>-1e-10);
  const bool hit_max (dl == max_len - r.length());

  // Calculate the column density. Because the density has a linear
  // gradient, this is equal to the density at the center of the
  // segment.
  
  // Note that we do the calculation of the gradient.ray in *AREPO
  // length units* because that requires the fewest operations

  assert(DP[dp_index()].task==ThisTask);

  if( (SphP[sphp_idx].Density>0) && 
      any(cell_->data()->get_absorber().densities()>0)) {

    // cell center of mass transformed to the dp_idx cell (in arepo units)
    const vec3d cell_center
      (DP[dp_index()].x-DP[primary_cells[dp_index()]].x+SphP[sphp_idx].Center[0],
       DP[dp_index()].y-DP[primary_cells[dp_index()]].y+SphP[sphp_idx].Center[1],
       DP[dp_index()].z-DP[primary_cells[dp_index()]].z+SphP[sphp_idx].Center[2]);

    const vec3d arepo_gradient(SphP[sphp_idx].Grad.drho[0],
			       SphP[sphp_idx].Grad.drho[1],
			       SphP[sphp_idx].Grad.drho[2]);

    vec3d pmid(r.position() + 0.5*dl*r.direction());
    pmid*=(1./lcon);

    const vec3d delta_center(pmid - cell_center);

    // This expression calculates the ratio between the density at the
    // segment midpoint to the mean density of the cell. Take care to
    // get all lengths in arepo units.
    T_float rho_multi 
      (1 + 
       dot(arepo_gradient, delta_center)/
       SphP[sphp_idx].Density);

#ifndef NDEBUG
    if(rho_multi<0) {
      // should never happen
      const int nearest1 = locate_brute(r.position()+1e-5*r.direction());
      const int nearest2 = locate_brute(r.position()+(dl-1e-5)*r.direction());
      const T_float rho1 = cell_density(dp_index(), r.position());
      const T_float rho2 = cell_density(dp_index(), r.position()+dl*r.direction());
      const T_float rho3 = cell_density(dp_index(), pmid);

      std::cerr << "Error: Rho_multi=" << rho_multi << " in cell " << dp_index() << " start point in " << nearest1 << " endpoint in " << nearest2 << std::endl;
      std::cerr << "Enterpoint rho=" << rho1 << " exitpoint rho=" << rho2 << " midpoint rho=" << rho3 << std::endl;
      std::cerr << r.position() << vec3d(r.position()+dl*r.direction()) << std::endl;
      std::cerr << "previous cell is " << previous_cell_ << std::endl;
      assert(rho_multi>=0);
    }
#endif

    // calculate the column density along the ray segment (can do this
    // in operator form but if only one density this is faster)
    if(dn.size()==1)
      dn(0) = cell_->data()->get_absorber().densities()(0) * rho_multi * dl;
    else
      dn = cell_->data()->get_absorber().densities() * rho_multi * dl;
    
    DEBUG(2,std::cout << "Determined dl=" << dl << ", rho_multi=" << rho_multi << ", cell density=" << cell_->data()->get_absorber().densities()(0) << ", column density=" << dn(0) << std::endl;);
  }
  else {
    DEBUG(2,std::cout << "Determined path length " << dl << " through empty cell" << std::endl;);
    // no need to do all this stuff if the cell density is zero.
    dn = 0;
  }

  
  if(hit_max) {
    // we hit max length, terminate tracker
    cell_=0;
    return dl;
  }

  // Update this cell tracker for next cell
  previous_cell_ = dp_index();
  dp_idx_ = next;
  if(next>=0)
    cell_ = &cell_->host_grid()->cell(sph_index(next));
  else
    cell_=0;

  DEBUG(2, if(cell_) std::cout << "Entering Voronoi " << dp_index() << std::endl;);
  DEBUG(1, if(!cell_) std::cout << "Exiting volume" << std::endl;);

  return dl;
}    

#endif
#endif
